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Abstract 


Gravimetric models of the geoid over Western Australia have been constructed using two 
adapted forms of Stokes’s integral; one uses the unmodified Stokes kernel and the other uses a 
deterministically modified kernel. These solutions use a combination of the complete expansion of 
the EGM96 global geopotential model with Australian gravity and terrain data. The resulting 
combined solutions for the geoid are compared with the control given by Global Positioning System 
(GPS) and Australian Height Datum heights at 63 points over Western Australia. The improved fit 
of the model that uses a modification to Stokes’s kernel indicates that this approach is more 
appropriate for gravimetric geoid computations over Western Australia. 


Introduction 


The geoid is the equipotential surface of the Earth’s 
gravity field, which corresponds most closely with mean 
sea-level (ignoring oceanographic effects) and undulates 
with respect to an oblate ellipsoidal model of the figure 
of the Earth. In 1849, GG Stokes published a solution to 
the geodetic boundary-value problem, which requires a 
global integration of surface gravity data over the Earth 
to compute the separation (N) between the geoid and 
reference ellipsoid (Stokes 1849). However, the 
incomplete global coverage and availability of accurate 
gravity measurements has precluded an exact 
determination of the geoid using Stokes’s formula. 
Instead, an approximate solution is used in practice, 
where only gravity data in and around the computation 
area are used. This approach is also attractive because of 
the increase in computational efficiency that is offered by 
working with a smaller integration area. 


In 1958, MS Molodensky (cited in Molodensky et al. 
1962) proposed a modification to Stokes’s formula to 
reduce the truncation error that results when gravity data 
are used over a limited area. However, Molodensky’s 
modification did not receive a great deal of attention at 
that time because of the contemporaneous availability of 
low-frequency global gravity field information, derived 
from the analysis of the orbits of artificial Earth satellites. 
These global geopotential models are expressed in terms 
of fully normalised spherical harmonic functions and are 
now routinely used in conjunction with terrestrial gravity 
data via a truncated form of Stokes’s integral (e.g. Vincent 
& Marsh 1973; Sideris & She 1995). This combined 
approach reduces the truncation error because its series 
expansion begins at a higher degree, where the truncation 
coefficients are smaller in magnitude (assuming that the 
global geopotential model is an exact fit to the low-degree 
terrestrial gravity field). Another advantage of this 
combined solution for the geoid is that it reduces the 
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impact of the spherical approximation inherent to the 
derivation of Stokes’s integral (e.g. Heiskanen & Moritz 
1967), the reason being that most of the geoid’s power is 
contained in the low frequencies. 


This paper will briefly review the approaches cur- 
rently used in regional gravimetric geoid computations 
and presents a compromise approach based on a high- 
degree global geopotential model and a low-degree 
modified Stokes kernel (Featherstone et al. 1998). An em- 
pirical comparison between gravimetric geoid models 
computed for Western Australia using the unmodified 
Stokes kernel and the compromise approach will be made 
with 63 discrete geoid heights. From these comparisons, 
it will be concluded that the compromise approach is 
more appropriate for gravimetric geoid computations in 
Western Australia. 


The generalised Stokes scheme 


A formal representation of the combination of a global 
geopotential model with terrestrial gravity data has been 
proposed by Vaníček & Sjoberg (1991), which they refer 
to as the generalised Stokes scheme for geoid computa- 
tion. Importantly, this satisfies a solution to the geodetic 
boundary-value problem when formulated for a higher 
than second-degree reference model of the figure of the 
Earth (Martinec & Vaníček, 1996). In this scheme, the low- 
frequency geoid undulations generated by a global 
geopotential model (N,,) are extended into the high fre- 
quencies by a global integration of complementary high- 
frequency, terrestrial gravity anomalies (Ag) using 


on n 
N = Ny + f SM (cosy) Ag™ siny dyda (1) 
0 0 


where k = R/4zy, R is the spherical Earth radius, y is 
normal gravity evaluated on the surface of the geocentric 
reference ellipsoid as required by Bruns’s formula (e.g. 
Heiskanen & Moritz 1967), y and a are the coordinates of 
spherical distance and azimuth angle about the 
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computation point, respectively, and S“(cos y) is the 
spheroidal form of Stokes’s integration kernel, which is 
implicit to the generalised scheme and has the series 
expansion 


S™ (cosy) = z 


n=M+1 


2n+1 
n-1 


P,(cosw) (2) 


where P (cos y) is the n degree Legendre polynomial. 


In Eq (1), the low-frequency component of the geoid 
(N,,) is computed from the fully normalised spherical 
harmonic coefficients that define the global geopotential 
model according to 


GM. a ar A xo —m 
3 F > (5) Y (6Cam cosmÀ + Sam sin mA) P? (cos 6) 
n=2 
(3) 


The corresponding high-frequency gravity anomalies 
(Ag™) are evaluated by subtracting the same spherical 
harmonic degrees of the same global geopotential model 
from the terrestrial gravity anomalies (Ag) according to 


CMe ( 


Nm 


m=0 


Ag! = Ag— oom 
(4) 
In Eqs (3) and (4), GM, is the product of the Newtonian 
gravitational constant and mass of the solid Earth, oceans 
and atmosphere, a is the equatorial radius of the reference 
ellipsoid, (r, 0, A) are the geocentric spherical_polar 
coordinates of each computation point, ÔC n and S,„„ are 
the fully normalised geopotential coefficients of degree n 
and order m, which have been reduced by the even zonal 
harmonics of the reference ellipsoid, and P „(cos 0) are 
the fully normalised associated Legendre functions. It is 
assumed that the reference ellipsoid is geocentric and has 
the same mass, potential and rotation rate as the geoid, 
such that the zero and first degree harmonics are 
inadmissible (e.g. Heiskanen & Moritz 1967). 


The degree (M) of reference spheroid chosen for the 
generalised Stokes scheme can be driven by the 
maximum degree of global geopotential model available, 
which is usually M „a = 360. However, there are more 
important considerations than simply taking the 
maximum degree of expansion available (e.g. 
Featherstone 1992). Firstly, the M_,, = 360 models are 
already combined solutions for the geoid because they 
are constructed from both satellite-derived and terrestrial 
gravity data. Therefore, the same terrestrial gravity data 
are usually used twice in Eq (1), which gives rise to 
unknown correlations between these data that are rarely 
accounted for or even acknowledged by some authors. 
Another consideration is the leakage of low-frequency 
errors from the terrestrial gravity anomalies into the 
combined solution for the geoid, which can be filtered by 
the spheroidal kernel due to the orthogonality of 
spherical harmonic functions over the sphere (e.g. 
Vaníček & Featherstone 1998). This is considered to be a 
desirable scenario because the low-frequency 
geopotential coefficients are currently the best source of 
this information, whereas terrestrial gravity data are 
subject to low-frequency errors. Therefore, choosing the 
degree of spheroid at, say, M = 20 in Eq (1), which is the 
limit of the reliable resolution of the satellite-derived 


n us ae -— 
5) (n-1) y (6C nm cos mÀ + Spm sin mA) P? (cos@) 
=? m=0 


geopotential coefficients, both avoids the correlations and 
reduces the leakage of terrestrial gravity anomaly errors. 


Reduction of the Approximation Error 


The generalised Stokes scheme 


The generalised Stokes scheme (Eq 1) remains subject 
to a truncation error when high-frequency terrestrial 
gravity anomalies are used over a limited area. 
Accordingly, there is an adjustment of Eq (1) that 
involves limiting the integration domain to a spherical 
cap of radius y (0 < y < 2) about each geoid computation 
point, which yields the approximation 


p 27 pvo 
N>=N; =Nu +n f l SM (cosy) Ag™ siny dy da 
o Jo 
(5) 


with a corresponding truncation error of 
2n pr 
6N, = sf i S™ (cosy) Ag™ siny dw da (6) 
0 vo 


such that N=N,+6N,. This truncation error can be 
expressed as a series expansion by 


ON, = 20K ye QM (to) Agn 


(7) 
n=M+1 
where the spheroidal truncation coefficients 
QM (vo) = | S™ (cosy) P, (cosy) siny dy (8) 
Wo 


can be evaluated using the algorithms of Paul (1973). The 
n degree surface spherical harmonic of the gravity 
anomaly can be evaluated for each n up to the maximum 
degree (M) of the geopotential model using 


GM say" om = 
Agee a (=) (n-1) 2 am cos mÀ + Spam sin mA) P% (cos 6) 
(9) 
As such, the truncation error in Eq (7) reduces to 
co 
Ni = 2r Y QM (0) Agn (10) 
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However, if Ag’ # 0 (2 < n < M), which is true if the 
global geopotential model is not an exact fit to the low- 
frequency terrestrial gravity anomalies, Eq (10) no longer 
holds and there is a leakage of any low-frequency gravity 
errors into the geoid solution when the integration is 
performed over a limited area (Vaníček & Featherstone 
1998). This is a direct consequence of the approximation 
of the generalised Stokes integral, which introduces the 
non-zero truncation coefficients in the region 2 <n <M 
because the orthogonality of spherical harmonics over the 
sphere breaks down under the approximation in Eq (5). 
Since Ag, only depend on the physical properties of the 
Earth, it becomes necessary to seek a modification to the 
integration kernel that reduces the magnitude of the 
truncation coefficients and hence the truncation error. 
Ideally, this modification should reduce the truncation 
error as well as adapting the kernel to behave as a 
(partial) high-pass filter and thus reduce the leakage of 
any low-frequency terrestrial gravity data errors into the 
geoid (Vaníček & Featherstone 1998). 


W E Featherstone: Gravimetric geoid models over Western Australia 


The remove-compute-restore scheme 


The remove-compute-restore technique (e.g. Torge 
1991) has almost become a standard approach in regional 
combined solutions for the geoid. However, most users 
of this approach make no attempt to modify the 
integration kernel and thus (further) reduce the 
truncation error or adapt its filtering properties. Instead, 
this scheme uses the unmodified kernel as originally 
introduced by Stokes (e.g. Sideris & She 1995; Smith & 
Milbert 1999). The remove-compute-restore approach 
also generally uses the maximum available degree of the 
global geopotential model (M a) which gives rise to the 
approximated data combination 


M 2r ppo 
N œ No = NMmas + cf S(cos) Ag™™= siny dy da 
o Jo 


(11) 
where the terms N „ma and Ag" are computed from the 
maximum available degree and order of a global 
geopotential model (Eqs. 3 and 4), and the unmodified 


Stokes kernel is given by 


co 


S(cos ý) = >» 


n=? 


2 1 
2E Pa(cost) 


= (12) 

In this combined solution for the geoid, the truncation 
error for degrees greater than M has to be neglected 
because it cannot be computed since the complete 
expansion of the global geopotential model has already 
been used, This also makes it subject to the correlation of 
errors between the global geopotential model and 
terrestrial gravity data used in the regional geoid 
solution. Moreover, no other attempt has been made to 
reduce the truncation error or adapt the filtering 
properties of the integration kernel. Admittedly, the 
truncation error will reduced a great deal if the global 
geopotential model is a good fit to the terrestrial gravity 
anomalies in the area of interest (i.e. if Ag” = 0 in the 
region 2 < n S Mpa): However, the penalty of taking this 
approach is that any errors in the terrestrial gravity 
anomalies to propagate unattenuated into the combined 
solution for the geoid (Vaníček & Featherstone 1998). lt 
must also be acknowledged that, despite these 
restrictions, the remove-compute-restore technique has 
delivered quite reasonable results (e.g. Sideris & She 1995; 
Smith & Milbert 1999). However, the question of whether 
a modified integration kernel in the generalised Stokes 
scheme will deliver even better results remains, and thus 
forms the primary aim of this investigation. 


Integration kernel modifications 


As argued above, it remains preferable to apply a 
modification to the approximated form of the generalised 
Stokes’s integral (Eq 5) and the remove-compute-restore 
approach (Eq 11) so as to further reduce the associated 
truncation error. Since Molodensky’s pioneering work, 
several other authors have proposed modifications to 
Stokes’s (1849) integral. These have been based on 
different criteria and can be broadly classified as 
deterministic modifications (e.g. Molodensky et al. 1962; 
Wong & Gore 1969; Meissl 1971; Heck & Griininger 1987; 
Vaníček & Kleusberg 1987; Vaníček & Sjöberg 1991; 
Featherstone ef al. 1998) and stochastic modifications (e.g. 
Wenzel 1982; Sjöberg 1991; Vaníček & Sjoberg 1991). The 
stochastic modifications, whilst offering an optimal 


L 
SM (cosy) = S™ (cosh) - $ (cos Uo) = y 2k 
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combination of the two data types together with a 
minimisation of the truncation error (in a least-squares 
sense), require reliable variance estimates of the data. 
However, the error characteristics of the terrestrial 
gravity data over Western Australia and over most other 
parts of the world are currently unknown, which renders 
the stochastic modifications of limited practical use. 
Therefore, the deterministic kernel modifications will 
have to be relied upon in the interim. 


The deterministic kernel modifications can be further 
divided into two broad categories; modifications that 
reduce the upper bound of the truncation error according 
to some prescribed norm, and modifications that improve 
the rate of convergence of the series expansion of the 
truncation error. The modification proposed by 
Featherstone et al. (1998) uses a combination of these, 
where the rate of convergence of the series expansion of 
an already-reduced truncation error by the L, norm 
(Vaníček & Kleusberg 1987) is accelerated from O(n?) to 
O(n?) through the approach proposed by Meissl (1971). 
This can be achieved either by setting the kernel to zero 
at the truncation radius through subtraction, or by 
choosing the truncation radius such that it coincides with 
a zero point of the kernel. 


The Featherstone ef al. (1998) modification is given by 


- tk(Vo) [Pk (cos Y) — P(cosyo)} 
(13) 


where the #,(y,) modification coefficients are computed 
from the solution of the following set of L-1 linear 
equations, once the truncation radius (y,) has been 
chosen 


k=2 


Pk 
aes 


k=2 


tk (Wo) enk (Vo) = QM (vo) (14) 


where the coefficients Q™ (y,) are given by Eq (8), and 


m 
€nk(Wo) = P,(cos Y) Py(cosw) sin ù dy (15) 
vo 

which can be evaluated using the recursive algorithms of 
Paul (1973). The degree of this kernel modification (L) 
can be chosen to be greater than, equal to or less than the 
degree of the reference spheroid (M) embedded in the 
generalised Stokes formula (Eq 1). If L > M, additional 
terms arise due to this disparate combination and should 
be computed or their omission acknowledged. 


A compromise of the combined solution 


The combined solution for the geoid considered in this 
study attempts to reach a compromise of the above two 
schemes, based on considerations of the data availability, 
their expected reliability, and a reduction of the trunca- 
tion error through the above deterministic modification 
of the generalised Stokes kernel. This compromise approach 
was used to compute the recent Australian gravimetric 
geoid model, AUSGeoid98 (Johnston & Featherstone 
1998). Mathematically, this is formalised as 


4 2n rpo 
N~N3= Ny... + nf | SM (cos) AgMn«= siny dy da 
o Jo 
(16) 
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where all terms are as defined earlier. Equation (16) 
utilises the maximum available expansion (M a) of the 
global geopotential model in conjunction with a low- 
degree (L) of deterministic kernel modification. This 
approach aims at reducing the truncation error so that it 
can be safely ignored, whilst relying more on the low- 
degree satellite solution by filtering a proportion of the 
low-frequency errors from the terrestrial gravity data (cf. 
Vaníček & Featherstone 1998). However, this choice is 
also driven by some practical considerations. Empirical 
studies by Featherstone (1992) indicate that the modified 
kernels become numerically unstable for large L and 
small y, which enforces a low degree of kernel 
modification when a small integration radius is used. For 
simplicity, the degree of kernel modification is chosen 
equal to the degree of spheroid used in the generalised 
scheme (i.e. L = M = 20). The integration radius was 
chosen to be y,=1°, since this value was empirically 
selected for AUSGeoid98 (Johnston & Featherstone 1998). 


It is argued that this compromise approach offers a 
geoid solution that is superior to the application of the 
remove-compute-restore technique with an unmodified 
kernel because of its further reduction of the truncation 
error and adaption of the filtering properties of the 
kernel. However, it is also important to acknowledge the 
deficiencies of this attempted compromise, which are the 
use of the high-frequencies in the global geopotential 
model (which can contain 80% noise; e.g. Lemoine et al. 
1998) and the correlations between the terrestrial gravity 
data in the region 20<n<M__. Therefore, empirical tests 
are used in Western Australia to determine whether the 
use of a deterministically modified integration kernel is 
more appropriate than using the unmodified kernel. 


Empirical tests in Western Australia 


The tests that follow compare four combined solutions 
for the geoid, computed using differing parameters in 
Eqs (11) and (16), with discrete geoid heights at 63 points 
across Western Australia. These control data are derived 
from Global Positioning System (GPS) measurements at 
Australian Height Datum (AHD) benchmarks. The 
difference between a GPS-derived ellipsoidal height and 
geodetically levelled height with respect to mean sea 
level gives a discrete estimate of the separation between 
the geoid and reference ellipsoid. Assuming that there 
are no errors in these contro] data, the gravimetric geoid 
solution that delivers the best fit can be considered to be 
the most suitable for Western Australia. 


Data and its processing 


The EGM96 global geopotential model (Lemoine et al. 
1998), complete to Mpa = 360, was used in this study. 
EGM9%6 is one of the most recent global geopotential 
models and was produced by the US National Imagery 
and Mapping Authority (NIMA), formerly the Defense 
Mapping Agency (DMA), and the US National 
Aeronautical and Space Administration’s (NASA) 
Goddard Space Flight Center (GSFC). Kirby et al. (1998) 
compared EGM96 with Australian gravity data 
(described below) and discrete geoid heights provided by 
co-located GPS and AHD data (also described below). 
This showed that EGM96 provides a slightly better 
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(though statistically not significant) representation of the 
Australian gravity field than its predecessors. 


The Australian Geological Survey Organisation’s 
(AGSO) 1996 national gravity data release, being the most 
recent available to the author, has been used in the 
computations. As most of the AGSO gravity data were 
collected and reduced predominantly for geophysical 
exploration purposes, they are not necessarily suited to 
the requirements of gravimetric geoid computation. As 
such, they have been validated according to the 
procedures in Featherstone et al. (1997) and the gravity 
anomalies computed using the more stringent geodetic 
approaches (e.g. Featherstone & Dentith 1997). 
Deficiencies were also found in the mean free-air gravity 
anomalies computed on land, due to the gravity data 
collection strategies used by AGSO. The land gravity 
observations were typically made along roads and tracks 
in areas of rugged terrain, which usually follow valleys 
or areas of least height variation. A similar, though 
opposite, effect occurs in central Australia where most 
observations were performed using helicopter transport. 
These observations were often located on raised ground 
where convenient landing spots were identified. 
Specifically, the computed mean gravity anomaly does 
not truly represent the actual mean gravity anomaly of 
an area. As such, the computed gravimetric geoid is 
biased by these gravity observation techniques. To 
numerically counter this biasing effect, a digital elevation 
model (described below), which gives a better 
representation of the mean terrain height than the gravity 
observation elevations, was used to reconstruct more 
representative mean gravity anomalies (Featherstone & 
Kirby 2000). 

Satellite-altimeter-derived marine gravity anomalies 
were used offshore Australia to supplement AGSO’s 
marine gravity data coverage. These gravity anomalies 
were computed from a combination of satellite-borne 
radar altimeter missions and supplied on a 2' by 2' 
geographical grid (Sandwell & Smith 1997). These data 
significantly improve the gravity data coverage offshore 
Australia and, as such, were expected to give an 
improved geoid solution in marine areas and on land 
areas close to the coast. However, errors thought to reside 
in the low-frequency altimeter-derived gravity anomalies 
adversely affect the gravimetric geoid near the coast, 
which was indicated by comparisons of the geoid 
solution with GPS and AHD data on land near the coast. 
The low-frequency errors in the satellite altimeter data 
are probably due to a combination of poorly-modelled 
near-shore sea-surface topography, tides and backscatter 
of altimeter returns from the land. Another explanation 
for this observed deficiency comes from the conversion 
of sea surface heights, which are the direct measurements 
taken by the altimeter, to gravity anomalies, then 
converting these gravity anomalies to geoid heights. At 
present, it is unknown how the errors, such as poorly- 
modelled sea-surface topography, propagate through 
these processes. Therefore, as an interim and practical 
solution, least squares collocation (e.g. Moritz 1980a) was 
used to ‘drape’ the altimeter gravity anomalies onto the 
AGSO marine gravity anomalies (Kirby & Forsberg 1998). 
This approach then improved the geoid solution on land 
near the coast, when compared with GPS and AHD data, 
over that achieved using only the AGSO marine gravity 


W E Featherstone: Gravimetric geoid models over Western Australia 


data or simple averaging of the AGSO marine gravity 
data and the satellite-altimeter data. 


Gravimetric terrain corrections, based on the national 
9" by 9" digital elevation model (Carrol & Morse 1996), 
were evaluated using Moritz’s (1968) formula via the fast 
Fourier transform or FFT (Kirby & Featherstone 1999). 
The FFT offers the only practical way to compute detailed 
terrain corrections on a continental scale, since prism 
integration at this scale could take several months to 
evaluate. lt was found by Kirby & Featherstone (1999) 
that the terrain correction computations had to be 
performed using a 27" by 27" grid to avoid the instability 
in Moritz’s terrain correction algorithm that occurs when 
using high-resolution digital terrain models close to the 
computation point (cf. Martinec et al. 1996). Associated 
with the gravimetric terrain correction are the primary 
and secondary indirect effects (e.g. Wichiencharoen 1982). 
The primary indirect effect accounts for the change in 
potential caused by the free-air gravity reduction and 
gravimetric terrain correction. The primary indirect effect 
was computed using the FFT on a 27" grid, which avoids 
the kernel instability and, moreover, is consistent with 
the terrain correction computations. The secondary 
indirect effect on gravity was computed by applying the 
free-air reduction over the geoid-compensated-geoid 
separation computed via the primary indirect effect. This 
resulted in an additional gravity term that was added to 
the gravity anomalies prior to geoid computation. 


A grid of residual gravity anomalies was computed 
from the AGSO gravity observations using the 
continuous curvature spline in tension algorithm (Smith 
& Wessel 1990; Wessel & Smith 1995). The term residual 
gravity anomalies is used to describe the terrain-corrected 
free-air gravity anomalies that have been reduced by the 
gravity anomalies implied by the complete M„„ = 360 
expansion of EGM96 (Eq 4). A regular grid of residual 
gravity anomalies is required for the residual geoid 
computations via the one-dimensional FFT technique. It 
is acknowledged that other gravity gridding techniques 
exist, such as least squares collocation (Moritz 1980a), but 
the continuous curvature spline in tension algorithm was 
readily available (Wessel & Smith 1995) and gives almost 
identical results in a considerably shorter computation 
time (Zhang 1997). For this study, a 2' by 2' grid of gravity 
anomalies was generated over the region by -11° < ¢< - 
37° and 110° < AS 131°, which covers the state of Western 
Australia. Table 1 summarises the statistical properties of 
the gravity anomalies and the residual gravity anomalies. 
The GRS80 reference ellipsoid (Moritz 1980b) has been 
used in all computations. Grey-scale images of the 
EGM96-implied gravity anomalies and residual gravity 
anomalies are shown in Figs 1 and 2, respectively. 


The statistical fit of the gravity anomalies implied by 
the EGM96 global geopotential model to the gravity field 
of Western Australia is poorer than that experienced in 
other parts of the world (cf. Forsberg & Featherstone 
1998). This is probably due to the larger uncertainty in 
the terrestrial gravity anomalies, which is caused princi- 
pally by errors in the gravity station elevations (e.g. 
Featherstone et al. 1997). However, since a large propor- 
tion of these data has been used in the construction of the 
EGM96 global geopotential model, a more likely explana- 
tion is an inaccuracy in the mean gravity anomalies used 
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Figure 1. The M a= 360 EGM96-implied gravity anomalies over 
Western Australia (units in mGal; Mercator’s projection from the 
GRS80 ellipsoid). 


112° 


-20° Fy 


24 y 


-28 9 


-32 f -32° 


36° 


112° : js 124° 128° 


Figure 2. The residual gravity anomalies over Western Australia 
(units in mGal; Mercator‘s projection from the GRS80 ellipsoid). 
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Table 1 


Minimum, maximum, mean, standard deviation (SD) and root mean square (RMS) for the gravity 


anomalies over Westem Australia (units in mGal). 


gravity anomalies Minimum Maximum Mean sD RMS 
terrain-corrected free-air 131.085 -147.678 -8.762 + 29.110 + 30.400 
residual 91.554 -101.236 0.009 + 10.968 + 10.968 


for this model. Lemoine et al. (1998) describe the con- 
struction of the JGP95E 5' by 5' digital elevation model, 
where the Australian 5' by 5' digital elevation model was 
used west of A = 140° and a National Imagery and Map- 
ping Agency (NIMA) digital elevation model was used 
east of A = 140°. The Australian 5' by 5’ digital elevation 
model used by NIMA was probably that constructed at 
the Australian National University from the elevations 
associated with the gravity station elevations. 
Featherstone & Kirby (2000) show that this digital eleva- 
tion model is biased because the gravity observations 
do not accurately represent the topographic morphology 
(see earlier discussion). Therefore, more reliance should 
be placed on the terrestrial gravity anomalies described 
earlier, hence the use of a low-degree kernel modifica- 
tion. 


Geoid computation by the 1D-FFT technique 


In the mid 1980s, the fast Fourier transform (FFT) 
technique began to find wide-spread use in gravimetric 
geoid computation, because of its efficient evaluation of 
convolution integrals when compared to quadrature- 
based numerical integration. For many years, the planar, 
two-dimensional FFT was used (e.g. Schwarz et al. 1990). 
Strang van Hees (1990) then introduced the spherical, 
two-dimensional FFT. However, both of these FFT 
approaches are subject to approximation errors, the most 
notable of which is the simplification of Stokes’s kernel. 
Forsberg & Sideris (1993) therefore proposed the 
spherical, multi-band FFT, which reduces the impact of 
the simplified kernel. Haagmans et al. (1993) then refined 
this approach to give the spherical, one-dimensional FFT, 
which requires no simplification of Stokes’s kernel. For 
this reason, the 1D-FFT has been used in this 
investigation so that the exact kernels in Eqs (12) and (13) 
can be used efficiently and without the need for any 
simplification. 

Another consideration is that remove-compute-restore 
determinations of the geoid using the 1D-FFT often 
convolve the whole rectangular grid of gravity anomalies 
over a region with the spherical Stokes kernel (e.g. Sideris 
& She 1995; Smith & Milbert 1999). On the other hand, 
quadrature-based geoid determinations only use gravity 
anomalies over a spherical integration radius about each 
computation point. Therefore, each approach results in a 
different truncation error due to the neglect of the 
residual gravity anomalies in the (different shaped) 
remote zones outside each integration domain. Both of 
these implementations are tested in this study, where in 
Eq (11) the spherical integration radius is set to y, = 7 so 
as to use the whole gravity data rectangle, and a limited 
spherical cap is used to mimic quadrature-based 
numerical integration. 


In order to make the 1D-FFT approach closely mimic 
quadrature-based numerical integration over a spherical 
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cap, the integration kernel is set to zero outside the 
truncation radius (y, = 1°) before transformation to the 
frequency domain. A further adaptation of this technique 
has been added to allow the 1D-FFT-based evaluation of 
Eq (16), where the modified kernel (Eq 13) was 
implemented by evaluating it before transformation to 
the frequency domain (Featherstone & Sideris 1998). 
Comparisons with quadrature-based numerical 
integration software that uses spherical caps and 
deterministically modified integration kernels 
(Featherstone, 1992) were used to verify these adaptations 
of the 1D-FFT. 


In what follows, the gravimetric geoid heights from 
each approach have been evaluated using the 1D-FFT for 


(i) the remove-compute-restore technique with an 
unmodified kernel and residual gravity anomalies over 
whole data area y = 7, 


(ii) the remove-compute-restore technique with an 
unmodified Stokes’s kernel and residual gravity 
anomalies over a limited spherical cap (y, = 1°), and 


(iii) the compromise approach in Eq (16) with a 
deterministically modified kernel and residual gravity 
anomalies over a limited spherical cap (y, = 1°). 


In the latter case, the degree of spheroid associated 
with the generalised Stokes scheme is M = 20 and the 
degree of deterministic kernel modification is the same. 
Of course, a large number of permutations of these 
parameters is possible by varying the degree of global 
geopotential model (M), integration radius (y,), degree of 
kernel modification (L), and even the type of kernel 
modification (see section Integration Kernel 
Modifications, above). However, only the above three 
cases are studied for Western Australia because the 
remove-compute-restore technique is used in many other 
parts of the world (e.g. Sideris & She 1995; Smith and 
Milbert 1999), and the L = M = 20 modification for a cap 
radius of y, = 1° was used for AUSGeoid98 (Johnston & 
Featherstone 1998). 


All geoid computations were conducted on a 2’ by 2' 
grid over an area bound by -12° < @<-36° and 112°<A< 
129°, which eliminates the edge effect associated with the 
one-degree integration radius (with the exception of the 
integration over the whole data area). It should be 
pointed out that this edge effect affects the whole 
computation area when the cap-radius is unlimited. 
Nevertheless, the geoid comparisons are conducted over 
the same area for the sake of consistency. 


Comparison of Geoid Results with GPS at AHD 
Benchmarks 


As is customary in almost all validations of 
gravimetric geoid models on land, the geoid results from 
this study were compared with Global Positioning 
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System (GPS) and geodetic levelling data to determine if 
any improvements are made when utilising a spherical 
integration cap and deterministically modified kernel. 
However, it must be pointed out this type of comparison 
is inevitably biased because the geodetic levelling data 
used in this investigation are based on the Australian 
Height Datum (AHD). The AHD is a normal orthometric 
height system based on estimates of mean sea level from 
30 tide gauges around Australia (Roelse et al. 1971). As 
such, it does not give an exact representation of the 
equipotential geoid (e.g. Featherstone 1998). Nevertheless, 
GPS and geodetic levelling data currently provide the 
only (partially) independent means with which to verify 
a gravimetric geoid model on land. Given that the 
primary geodetic application of a gravimetric geoid 
model is to transform GPS-derived ellipsoidal heights to 
the AHD, this type of analysis is also useful to ascertain 
the gravimetric geoid models’ performance for this 
application. 


For each of the three cases investigated, the 
gravimetric geoid solution was bi-cubically interpolated 
from the 2' by 2' grid and statistically compared with 63 
discrete geoid heights. These were derived geometrically 
from the most precise GPS networks available in Western 
Australia (Morgan et al. 1996; Stewart et al. 1997) at points 
that are co-located with geodetically levelled heights of 
third-order, or better, on the AHD. Table 2 shows a 
statistical summary of the differences between the 63 
control geoid heights and the results from the M_,=360 
expansion of EGM96 alone (Eq 3), the 1D-FFT 
implementations of Eq (11) with y = zand y = 1°, and 
Eq (16) with y, = 1°. The mean and root mean square 
(RMS) differences in Table 2 should not be relied upon to 
choose the most accurate geoid solution because any 
gravimetric determination of the geoid is deficient in 
scale (i.e. only the gravimetric geoid’s precision can be 
estimated). This scale deficiency due to the imperfect 
knowledge of the mass of the Earth [cf the text 
immediately after Eqs (3) and (4)]. Accordingly, only the 
standard deviations of the fit of each gravimetric geoid 
model to the control data should be used to assess the 
performance in terms of precision of each model; the 
mean and RMS values are only included for the sake of 
completeness. 


The difference between the fit of each geoid model in 
Table 2 is not always significant in a statistical sense, 
when considering that the random error budget of the 
GPS data is ~0.05 m (Morgan et al. 1996; Stewart ef al. 
1997) and distortions of the AHD from the geoid are of 
the order of 1 m (Roelse et al. 1971; Featherstone & 
Stewart 1998; Featherstone 1998). Despite these caveats, 


some useful inferences can be made from these results as 
follows. 


Firstly, the use of the whole gravity data area in the 
combined solution for the geoid (Eq 11 with y, = z) gives 
a result that is worse than using the EGM96 global 
geopotential model alone (Table 2). Whilst the use of the 
whole data area appears to be appropriate in other parts 
of the world (e.g. Sideris & She 1995; Smith & Milbert 
1999), a considerably worse result is achieved in Australia 
when using this approach (Table 2; Forsberg & 
Featherstone 1998; Johnston & Featherstone 1998). This is 
a slightly unexpected result, since including more data in 
the geoid solution should yield a better result. This is 
either due to noise in the Australian gravity anomalies, 
which been estimated to be approximately + 2 mGal 
(Featherstone et al. 1997), or the theoretical basis of this 
approach is unsuitable for Australia. 


Next, from Table 2 there is an improvement over the 
results of using only the EGM96 global geopotential 


m2 116° 120° 124° 128° 


Figure 3. The gravimetric geoid over Western Australia, com- 
puted using Eq (16) with y, = 1° and S,,""(cos y); units in metres; 
Mercator’s projection from the GRS80 ellipsoid. 


Table 2 


Minimum, maximum, mean, standard deviation (SD) and root mean square (RMS) for the absolute 
between the 63 control GPS-AHD heights and the gravimetric geoid heights computed from EGM96 
only, Eq (11) with y = wand y = 1° and Eq (16) with y = 1° (units in metres). 


Minimum 
M „a = 360 expansion of EGM96 (Eq 3) 0.889 
Eq (11) with y = p and S(cos y) 0.751 
Eq (11) with y = 1° and S(cos y) 0.625 
Eq (11) with y = 1° and S,,*(cos y) 0.761 
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Maximum Mean SD RMS 
-0.328 0.148 +0.274 +0.311 
-0.776 -0.092 T 0:3620 R0374 
-0.664 0.118 E0249 +0.276 
-0.308 0.140 EEO meet 0263 
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120° 


112° 128° 


Figure 4. The differences between the 63 GPS-AHD discrete 
geoid heights (circles) and the gravimetric geoid (Fig 3); units in 
metres; Mercator’s projection from the GRS80 ellipsoid. 


model offered by the use of a limited integration radius 
in Eq (11). Importantly, there is a further improvement 
offered by using the deterministically modified 
integration kernel in Eq (16) for the same integration 
radius. Though inconclusive due to the perceived 
uncertainties in the control data, these results do indicate 
that the use of deterministically modified integration 
kernel is a preferable approach for geoid computations in 
Western Australia. This is bearing in mind that the 
principal application of the gravimetric geoid in Australia 
is to transform GPS-derived ellipsoidal heights to the 
AHD. Fig 3 shows a contour map of the gravimetric 
geoid model of Western Australia that has been 
computed using Eq (16) with y, = 1° and L = M = 20. 


Finally, it is informative to view the geographical dis- 
tribution of the discrepancies between the 63 GPS-AHD 
discrete geoid heights and the gravimetric geoid com- 
puted from the modified Stokes integral (Fig 4). There are 
systematic differences between the gravimetric geoid and 
the control data, where these features (highs, lows and 
slopes) are defined by more than one control point. Given 
the relatively small error budget of the GPS data (= 0.05 
m), it becomes difficult to ascertain whether these sys- 
tematic differences are due to distortions in the AHD, the 
gravimetric geoid, or both. It is currently impossible to 
isolate the exact source, but following the arguments in 
Featherstone & Stewart (1998), it is more likely that these 
errors are due to distortions in the AHD. These are intro- 
duced because the AHD is constrained to mean sea level 
observed over two years at 30 tide gauges around Aus- 
tralia (Roelse et al. 1971) and thus does not necessarily 
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represent an equipotential surface of the Earth’s gravity 
field (Featherstone 1998). However, the possibility of sys- 
tematic errors in the gravimetric geoid cannot be ruled 
out as an explanation for these discrepancies. 


Conclusion 


From these gravimetric geoid results over Western 
Australia, it is clear that 1D-FFT geoid computations 
should use a spherical cap of limited extent instead of the 
whole gravity data grid. Small (though not statistically 
significant when considering the errors in the GPS and 
AHD data) improvements are also observed when a 
deterministically modified integration kernel is used over 
a spherical cap in the compromise approach to 
gravimetric geoid determination (Eq 16). This 
improvement is probably because the kernel modification 
reduces the truncation error so that its neglect has less 
impact on the solution and adapts the filtering properties 
of the kernel thereby reducing the effect of low-frequency 
terrestrial gravity anomaly errors. Therefore, it is 
recommended that limited integration caps and modified 
Stokes’s kernels are used in spectral geoid determinations 
of Western Australia. 
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